library(tidyverse)
library(ggpubr)
library(plotly)

upper_to_int <- function(x) map_vec(x, utf8ToInt) - 64

upper_to_int('A')

# load data ----
taqman <- read_tsv('00_util_scripts/data/240821-taqman-g2r.txt', skip = 1) |>
  mutate(ratio = `465-510`/`533-580`)

taqman <- taqman |>
  filter(Call != 'Negative')

taqman |>
  count(Call)

taqman <- taqman |>
  mutate(col = str_remove(Pos, '[A-Z]') |> as.numeric(),
         row = str_extract(Pos, '[A-Z]') |> upper_to_int(),
         id = 1 + 8 * (col) - row / 2) |>
  filter(id <= 77)

# roche calculation result
taqman |>
  ggplot(aes(`465-510`, `533-580`, color = Call)) +
  geom_point()

taqman <- taqman |>
  mutate(genotype = case_match(Call,
                               'Both Alleles' ~ 'GG',
                               'Allele Y' ~ 'GR',
                               c('Unknown','') ~ 'RR',
                               .default = 'Negative control') |>
           fct_relevel('Negative control'),
         sample = ifelse(id <= 10, 'PBMC', 'normal')) |>
  arrange(id) |>
  write_tsv('taqman0822.i2t.tsv')
  
taqman |>
  ggplot(aes(`465-510`, `533-580`, color = genotype, shape = sample))+
  geom_point(size = 3)+
  theme_pubr() +
  labs_pubr() +
  expand_limits(x = 0, y = 0) +
  labs(title = 'Taqman result of IGHG1-G396R genotype from 77 samples')

ggplotly()

unknown <- taqman |>
  filter(Pos %in% c('B16','F19','C11'))

taqman <- taqman |>
  filter(!(Pos %in% c('B16','F19','C11')))

taqman |>
  ggplot(aes(row, col, fill = Call, label = Name)) +
  geom_raster() +
  geom_text()

# kmeans results are randomly produced ------
# so set nstart to run multiple times and select the best clustering
kmns_res <- taqman |>
  filter(Call != 'Negative', Pos != 'J11') |>
  select(ratio) %>%
  kmeans(centers = 3, nstart = 25) %>%
  pluck("cluster") %>%
  as.character()

table(kmns_res)

taqman_res <- taqman %>%
  filter(Call != 'Negative', Pos != 'J11') |>
  mutate(kmeans = kmns_res)

taqman_res |>
  ggplot(aes("",ratio, color = kmeans)) +
  geom_point()

taqman_res <- taqman_res %>%
  group_by(kmeans) %>%
  summarise(order = mean(ratio)) %>%
  arrange(order) %>%
  add_column(genotype = c("II", "IT", "TT")) %>%
  select(-order) %>%
  right_join(taqman_res, multiple = 'all')

unknown$genotype <- 'unknown'

taqman_res <- unknown |> bind_rows(taqman_res)

taqman_res %>%
  ggplot(aes(`465-510`, `533-580`, text = Pos))+
  geom_point(aes(color = genotype), size = 3)+
  theme_pubr() +
  labs_pubr() +
  expand_limits(x = 0, y = 0)

ggplotly()

write_csv(taqman_res, 'mx-0609-taqman.csv')

taqman_res |>
  mutate(Name = 1:39) |>
  ggplot(aes(row, col, fill = genotype, label = Name)) +
  geom_raster() +
  geom_text()

taqman_res <- taqman_res |>
  mutate(kmeans = case_match(genotype, 'GG' ~ 1, 'GR' ~ 2, .default = 3))

write_csv(taqman_res, 'tmp_taqman.csv')

taqman_old <- read_csv('tmp_taqman.csv')

# merge two results -----
taqman_merge <- taqman_old |>
  left_join(taqman_res,by = join_by(Name, Pos, row, col), suffix = c(".old",".new"))

taqman_merge |>
  ggplot(aes(ratio.old, ratio.new, size = Score.old + Score.new)) +
  geom_point() +
  stat_cor()

write_csv(taqman_merge, 'taqman_merge.csv')

taqman_merge %>%
  ggplot(aes(x = genotype.new))+
  geom_bar(aes(fill = genotype.old))+
  geom_text(aes(label = after_stat(count)), stat='count', vjust = 1.5, colour = "white", size = 7)+
  theme_pubr() +
  labs_pubr() +
  theme(text = element_text(size=16))

taqman_merge %>%
  ggplot(aes(x = Call.new))+
  geom_bar(aes(fill = Call.old))+
  geom_text(aes(label = after_stat(count)), stat='count', vjust = 1.5, colour = "white", size = 7)+
  theme_pubr() +
  labs_pubr() +
  theme(text = element_text(size=16))

taqman_additive <- taqman_merge |>
  mutate(vic = `465-510.old` + `465-510.new`, fam = `533-580.old` + `533-580.new`, ratio.add = vic/fam)

kmns_res <- taqman_additive %>%
  select(ratio.add) %>%
  kmeans(centers = 3, nstart = 25) %>%
  pluck("cluster") %>%
  as.character()

table(kmns_res)

taqman_additive <- taqman_additive |>
  mutate(kmeans.add = kmns_res)

taqman_additive |>
  ggplot(aes(x = '', y = ratio.add, color = kmeans.add)) +
  geom_point()

write_csv(taqman_additive, 'taqman_additive.csv')

taqman_additive |>
  ggplot(aes(vic, fam, color = kmeans.add, text = Name)) +
  geom_point() +
  scale_color_discrete(label = c('GG','GR','RR')) +
  labs(x = 'VIC', y = 'FAM', color = 'IGHG1-G396R genotype') +
  theme_pubr() +
  labs_pubr()

viz <- last_plot()

# interactive plot 
plotly::ggplotly(viz, tooltip = c("text"))

taqman_additive %>%
  mutate(kmeans.add = as.numeric(kmeans.add) |> case_match(1 ~ 'GG', 2 ~ 'GR', .default = 'RR')) |>
  ggplot(aes(x = kmeans.add))+
  geom_bar(aes(fill = kmeans.add))+
  geom_text(aes(label = after_stat(count)), stat='count', vjust = 1.5, colour = "white", size = 7)+
  theme_pubr() +
  labs_pubr() +
  labs(x = 'IGHG1-G396R genotype') +
  theme(text = element_text(size=16), legend.position = 'none')

viz2 <- last_plot()
viz + viz2 + patchwork::plot_layout(guides = 'keep')
